#!/bin/bash
#
echo ''
echo "`basename $0` ver 1.0"
echo 'Script to fit the EOS from a VASP run.'
echo "Usage:"
echo "   `basename $0` [eng_a.dat]"
echo 'The first column should be volume/atom in unit of Ang^3/atom;'
echo 'the second column should be energy/atom in unit of eV/atom.'
echo ''
echo 'Author: LT Kong, konglt@sjtu.edu.cn'
echo ''
#
evf="eng_a.dat"
if [ "$#" -ge 1 ]; then
   evf=$1
fi
if [ ! -f "${evf}" ]; then
   echo "File: ${evf} not found!"
   exit 2
fi
#
mineng=`sort -n -k 2 ${evf}|head -1`
E0=`echo $mineng|awk '{print $2}'`
V0=`echo $mineng|awk '{print $1}'`
#
cat > .vasp_ev_fit_scpt << EOF

EOS(x) = B0*x/Bp * ( (V0/x)**(Bp)/(Bp-1.) + 1.) - B0*V0/(Bp-1.) + E0
#
set fit errorvariables
#
B0 = 1.; Bp = 4.; E0 = $E0; V0 = $V0
fit EOS(x) "${evf}" u (\$1):(\$2) via B0, Bp, V0, E0

set term post enha colo 20
set out 'eos_fit.eps'

set xlabel 'volume'
set ylabel 'energy'

plot EOS(x) w l lw 2 t 'EOS', "${evf}" u (\$1):(\$2) w p pt 7 t 'VASP'

unset out
set term x11
print ''
print 'B0 = ', B0*160.217733,' +/- ',B0_err*160.217733, ' GPa'
print 'E0 = ', E0,' +/- ', E0_err, ' eV/atom'
print 'Bp = ', Bp, ' +/- ', Bp_err
print 'V0 = ', V0, ' +/- ', V0_err, ' Ang^3/atom'
print ''
print 'If FCC, a_0 = ',(V0*4.)**(1./3.), ' Ang'
print 'If BCC, a_0 = ',(V0*2.)**(1./3.), ' Ang'
print ''

# output a png file
set term png enha trans truec
set out "eos_fit.png"
replot

EOF
#
gnuplot .vasp_ev_fit_scpt
evince eos_fit.eps
rm -rf .vasp_ev_fit_scpt fit.log

exit 0
